#=========================================#
# Replication file for 
# High-level Visit and National Security Policy: 
# Evidence from A Quasi-Experiment in Taiwan
# Nov 9 2022
# austin.wang@unlv.edu

#clean the environment
rm(list=ls(all=TRUE))

#summarySE is a useful function for descriptive analysis
summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
  library(plyr)
  length2 <- function (x, na.rm=FALSE) {
    if (na.rm) sum(!is.na(x))
    else       length(x)
  }
  datac <- ddply(data, groupvars, .drop=.drop,
                 .fun = function(xx, col) {
                   c(N    = length2(xx[[col]], na.rm=na.rm),
                     mean = mean   (xx[[col]], na.rm=na.rm),
                     sd   = sd     (xx[[col]], na.rm=na.rm)
                   )
                 },
                 measurevar
  )
  
  datac <- rename(datac, c("mean" = measurevar))
  
  datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean
  
  ciMult <- qt(conf.interval/2 + .5, datac$N-1)
  datac$ci <- datac$se * ciMult
  
  return(datac)
}


#Use left, mid, and right for creating the time variable.
left = function(text, num_char) {
  substr(text, 1, num_char)
}

mid = function(text, start_num, num_char) {
  substr(text, start_num, start_num + num_char - 1)
}

right = function(text, num_char) {
  substr(text, nchar(text) - (num_char-1), nchar(text))
}


#Calling library
library(foreign)
library(car)              
library(ggplot2)

#read the data from the same folder, if needed
fullpath=rstudioapi::getSourceEditorContext()$path   #Make sure where is your working file location
currentfilepath=substr(fullpath,1,nchar(fullpath)-35) #35=length of file name (including.r)
currentfilepath
setwd(currentfilepath)
project1 = read.csv("PP2197E56.csv")
head(project1)
summary(project1)


#Please revise the link to your local address
#project1 = read.csv("C:/Users/weary/Dropbox/Chen_Wang_Wu_Yeh_Projects/National Security Surveys/High-level Visit and War Support/PP2197E56.csv",header=T,fileEncoding="UTF-8-BOM")

#Project 1 should have 1500 obs and 120 variables.

#Month variable
project1$month=left(project1$startdate,1)
table(project1$month)
#51 in May and 1449 in June

#Date variable
project1$date=NA

for(i in 1:1500){
  if(project1$month[i]==6){
    project1$date[i]=right(left(project1$startdate[i],3),1)
  }
  if(project1$month[i]==5){
    project1$date[i]=right(left(project1$startdate[i],4),2)
  }
}
project1$day=left(project1$startdate,1)

table(project1$date)
project1$date[project1$date==31]=0
project1$date=as.numeric(project1$date)
table(project1$date)

#51 in May 31
#660 in June 1...
#and 237 in June 8.

#Create the senator visit variable Davos
#Visited in the morning of June 6, so binary coded.
project1$Davos=0
project1$Davos[project1$date>=6]=1
table(project1$Davos)
#400 surveyed after the visit, 1100 before the visit 
#P.8-p.9 in the Data Collection Paragraph



#Control variables, Appendix Table A1
#Taiwan Identity, binary coded, Appendix Table A1
project1$TWIdentity=recode(project1$Q28,"1=1;else=0")
table(project1$TWIdentity)
#gender, Appendix Table A1
project1$female=NA
project1$female=recode(project1$Q57,"1=0;2=1")
table(project1$female)
#age, Appendix Table A1
project1$age=NA
project1$age=110-project1$Q29
table(project1$age)
#education, Appendix Table A1
project1$edu=project1$Q30
table(project1$edu)
#mainlander, Appendix Table A1
project1$mainlander=0
project1$mainlander[project1$Q31==3]=1
table(project1$mainlander)

#Party Identity, Appendix Table A1
project1$blue=0
project1$blue[project1$partyid2==1]=1
project1$blue[project1$partyid2==3]=1
project1$blue[project1$partyid2==5]=1
project1$blue[project1$partyid2==37]=1

project1$green=0
project1$green[project1$partyid2==2]=1
project1$green[project1$partyid2==6]=1
project1$green[project1$partyid2==20]=1
project1$green[project1$partyid2==21]=1
project1$green[project1$partyid2==23]=1
project1$green[project1$partyid2==36]=1
table(project1$blue)
table(project1$green)


#t.test and chisquare test, Appendix Table A2
t.test(age~Davos,data=project1)   #sig. 
t.test(edu~Davos,data=project1)   #sig.
chisq.test(table(project1$TWIdentity,project1$Davos))  #sig.
chisq.test(table(project1$female,project1$Davos))  #not sig
chisq.test(table(project1$blue,project1$Davos))  #sig
chisq.test(table(project1$green,project1$Davos))  #not sig

#==================================================#
#            Dependent variable                    #
#==================================================#


#Data collection, p.8

#Strong
#Do you agree with the statement that ��Only a stronger national defense can maintain peace"
project1$Enhance=recode(project1$Q15,"1=2;2=1;3=-1;4=-2")
table(project1$Q15,project1$Enhance)
summary(project1$Enhance)

#Confidence
#Do you think that our military is stronger than before
project1$strong=recode(project1$Q19,"1=1;else=0")
table(project1$Q19,project1$strong)

#Approve
#Are you satisfied with the current government��s overall performance in national defense
project1$satisfy=recode(project1$Q22,"1=2;2=1;3=-1;4=-2")
table(project1$Q22,project1$satisfy)


#t.test without control, Appendix Table A4
t.test(project1$Enhance[project1$Davos==0],
       project1$Enhance[project1$Davos==1])  #p=0.067
t.test(project1$strong[project1$Davos==0],
       project1$strong[project1$Davos==1])  #p<0.001
t.test(project1$satisfy[project1$Davos==0],
       project1$satisfy[project1$Davos==1])  #p<0.001


#=================================================#
#         Propensity Score Matching, p.10         #
#=================================================#
library(MatchIt)
match_output=matchit(Davos~blue+green+edu+female+age+TWIdentity,data=project1,
                     method="nearest",ratio=1,replace=FALSE)
summary(match_output)  #All 400 in the treatment group is paired
plot(match_output,type="jitter")  #Appendix Figure A3
plot(match_output,type="hist")  #Appendix Figure A3

plot(summary(match_output),xlim=c(0,2))  #Appendix Figure A3 - very good matched result

#create the output file after matching
m.data1 = match.data(match_output) 

##In the end, we can compare the two groups directly

#==============Table 1, P.10===============#
#Strong
t.test(m.data1$Enhance[m.data1$Davos==0],
       m.data1$Enhance[m.data1$Davos==1]) #p=0.046
#Confidence 
t.test(m.data1$strong[m.data1$Davos==0],
       m.data1$strong[m.data1$Davos==1])   #p=0.02 around, because of matching

#satisfy  
t.test(m.data1$satisfy[m.data1$Davos==0],
       m.data1$satisfy[m.data1$Davos==1])  #p=0.0609


#==============Figure 1, P.11===============#
kk1=data.frame(summarySE(m.data1,"Enhance","Davos"))
kk2=data.frame(summarySE(m.data1,"strong","Davos"))
kk3=data.frame(summarySE(m.data1,"satisfy","Davos"))

#combine the three summaries together
kk1[c(3:6),]=NA
kk1[3:4,]=kk2
kk1[5:6,]=kk3

kk1$group2=c("Before visit","After visit","Before visit","After visit","Before visit","After visit")
kk1$var=c("V1: Strong","V1: Strong",
          "V2: Support","V2: Support",
          "V3: Approval","V3: Approval")
kk1$group2=relevel(as.factor(kk1$group2), ref="Before visit")
library(ggsignif) 
ggplot(data=kk1,aes(x=var,y=Enhance, fill = group2))+geom_col(position = "dodge")+theme_bw()+
  geom_errorbar(aes(ymin=Enhance-ci, ymax=Enhance+ci), width=.2,position = position_dodge(width=0.9))+
  ylab("Mean value of Dependent Variables")+xlab("")+scale_fill_manual(values=c("#999999", "#56B4E9"))+
  theme(legend.title = element_blank(),legend.position = c(0.9, 0.9))+ylim(-0.5,2.5)+
  geom_signif(
    y_position = c(1.75, 1.75, 1.75), xmin = c(0.8, 1.8, 2.8), xmax = c(1.2, 2.2, 3.2),
    annotation = c("p=0.046*", "p=0.023*", "p=0.061+"), tip_length = 0.005, textsize = 4, size = 0.7 )  



#=================================================#
#         Regression Continuity Design, p.11      #
#=================================================#


#=================Table 2, p.12=======================#

model1=lm(Enhance~Davos,data=project1)
summary(model1)
model2=lm(Enhance~Davos+female+age+edu+blue+green+TWIdentity,data=project1)

model3=glm(strong~Davos,data=project1, 
           family = binomial(link = "logit"))
model4=glm(strong~Davos+female+age+edu+blue+green+TWIdentity,data=project1, 
           family = binomial(link = "logit"))

model3=glm(strong~Davos,data=project1, 
           family = binomial(link = "logit"))
model4=glm(strong~Davos+female+age+edu+blue+green+TWIdentity,data=project1, 
           family = binomial(link = "logit"))

model5=lm(satisfy~Davos,data=project1)
model6=lm(satisfy~Davos+female+age+edu+blue+green+TWIdentity,data=project1)

library(stargazer)

stargazer(model1,model2,model3,model4,model5,model6,type="text")


#============Table 3, p.14, no heterogeneity================#

model11=lm(Enhance~Davos+
             Davos*blue+Davos*green+TWIdentity+edu+female+age,data=project1)
model12=lm(Enhance~Davos+
             blue+green+Davos*TWIdentity+edu+female+age,data=project1)
model21=glm(strong~Davos+
              Davos*blue+Davos*green+TWIdentity+edu+female+age,data=project1, 
            family = binomial(link = "logit"))
model22=glm(strong~Davos+
              blue+green+ Davos*TWIdentity+edu+female+age,data=project1, 
            family = binomial(link = "logit"))
model31=lm(satisfy~Davos+
             Davos*blue+Davos*green+TWIdentity+edu+female+age,data=project1)
model32=lm(satisfy~Davos+
             blue+green+Davos*TWIdentity+edu+female+age,data=project1)
stargazer(model11,model12, model21,model22, model31,model32,type="text")


#============Appendix Figure 5, p.12================#
model2=lm(Enhance~as.factor(date)+
            blue+green+edu+female+age+TWIdentity,data=project1)
summary(model2)
library(sjPlot)
plot_model(model2,type="pred",term="date")+theme_bw()+
  geom_line()+
  geom_vline(xintercept=5.5,lty=2)+
  annotate(geom="text", x=5.5, y=0.9, label="Visit at 8am, June 6")+
  ylab("Agree that Enhancing military power to maintain peace (-2 - +2)")+xlab("Survey date (May 31-June 8, 2021)")+  ggtitle("")

model4=glm(strong~as.factor(date)+
             blue+green+edu+female+age+TWIdentity,data=project1, 
           family = binomial(link = "logit"))
plot_model(model4,type="pred",term="date")+theme_bw()+
  geom_line()+
  geom_vline(xintercept=5.5,lty=2)+
  annotate(geom="text", x=5.5, y=0.15, label="Visit at 8am, June 6")+
  ylab("% Believe that the military is stronger than before (0/1)")+xlab("Survey date (May 31-June 8, 2021)")+
  ggtitle("")

model6=lm(satisfy~as.factor(date)+
            blue+green+edu+female+age+TWIdentity,data=project1)
plot_model(model6,type="pred",term="date")+theme_bw()+
  geom_line()+
  geom_vline(xintercept=5.5,lty=2)+
  annotate(geom="text", x=5.5, y=-0.5, label="Visit at 8am, June 6")+
  ylab("Satisfy Defense Performance (-2 - +2)")+xlab("Survey date (May 31-June 8, 2021)")+
  ggtitle("")
